clear all 

folderfig = '/Users/wittwer/Dropbox/5_Jason and Milena/data/Matlab_JMP_2R1/main_estimation_submission2R1/3_figures/';
foldermats= '/Users/wittwer/Dropbox/5_Jason and Milena/data/Matlab_JMP_2R1/main_estimation_submission2R1/2_output/';


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%CHOOSE SPECIFICATION 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

zerobenefit =0; %--> CAREFUL: most parameters are automatically chosen below
    %=1 only if using the homedealer model but shutting off the loyalty benefit - benchmark model with fixed base qualitites
    %=0 else

homedealer=1;
eta=0; 
quantityrobus=0;
keepall=1;
day=1;
fej=1;
loyalty=1;
largetrade=0;tradesize=0; % or largetrade=1; tradesize=0;

if zerobenefit ==1 
    homedealer=0; day=1; keepall=1; largetrade=0; tradesize=0; quantityrobusL=0;
end 

if homedealer==1
        xij_for_investor=0;  
    else 
        xij_for_investor=1;  
end

%fix quotes as observed in the data
flexible=0;
    %=1, quotes as in equilibrium
    %=0, quotes as observed in the data

sideshow=1;
if sideshow==1
    disp(['BUY-SIDE'])
else 
    disp(['SELL-SIDE'])
end 

retailer=0;
side=sideshow;
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%CREATE MODEL FIT FIGURES - FIGURE 9 and APPENDIX FIGURE A6
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if zerobenefit==0
        temp=['/counterfac_SQ_',[num2str(side),'_' num2str(retailer), '_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(tradesize), '_', num2str(fej), '_', num2str(loyalty)]]; 
 else
        temp=['/counterfac_SQ_',[num2str(side),'_' num2str(retailer), '_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(tradesize), '_', num2str(fej), '_', num2str(loyalty), '_zerobenefit']]; 
        homedealer=0;
end
load(fullfile([foldermats,temp,'.mat']))


disp(['TRADING PROBABILITIES--------------------------------------------------------------------------'])

%Model
input = nanmean(sHjSQ + rhojSQ,2)*100;
disp(['MODEL: Average probability to trade with the home dealer (on or off the platform) is: ', num2str(nanmean(nanmean(sHjSQ + rhojSQ,2))*100), ' in CI = ', num2str(CI_fun(input))])
plotinput1a = input; 

%Data:
data = Y(Y(:,4)==retailer & Y(:,5)==side,:); 
share_H = zeros(T,max(Dmax));
for t=1:T
    for j=1:max(Dmax)          
            N_H = sum(data(data(:,3)==t & data(:,17)==j,26));
            N_T = sum(data(:,3)==t & data(:,17)==j);
            share_H(t,j) = N_H/N_T;
    end 
end 
input = nanmean(share_H,2)*100;
disp(['DATA: Average probability to trade with the home dealer (on or off the platform) is: ', num2str(nanmean(input)), ' in CI = ', num2str(CI_fun(input))])
plotinput1b = input; 
disp(['----------'])  

%Model
input = 100;
disp(['MODEL: Average probability to trade with the home dealer (off the platform) is: ', num2str(100), ' in CI = ', num2str(CI_fun(input))])

%Data:
data = Y(Y(:,4)==retailer & Y(:,5)==side,:); 
rhojSQ2_data = zeros(T,max(Dmax));
for t=1:T
    for j=1:max(Dmax)          
            N_H = sum(data(data(:,3)==t & data(:,17)==j & data(:,10)==1,26));
            N_T = sum(data(:,3)==t & data(:,17)==j & data(:,10)==1);
            rhojSQ2_data(t,j) = N_H/N_T;
    end 
end  
input = nanmean(rhojSQ2_data,2)*100;
disp(['DATA: Average probability to trade with the home dealer (off the platform) is: ', num2str(nanmean(input)), ' in CI = ', num2str(CI_fun(input))])
disp(['----------'])      
   

%Model:
disp(['MODEL: Average probability for R to trade with the home dealer (off the platform) is: ', num2str(100), ' in CI = ', num2str(CI_fun(100))])
plotinput2a = 100; 

%Data:
data = Y(Y(:,4)==1 & Y(:,5)==side,:); 
rhojSQ2_data = zeros(T,max(Dmax));
for t=1:T
    for j=1:max(Dmax)          
            N_H = sum(data(data(:,3)==t & data(:,17)==j & data(:,10)==1,26));
            N_T = sum(data(:,3)==t & data(:,17)==j & data(:,10)==1);
            rhojSQ2_data(t,j) = N_H/N_T;
    end 
end  
input = nanmean(rhojSQ2_data,2)*100;
disp(['DATA: Average probability for R to trade with the home dealer (off the platform) is: ', num2str(nanmean(input)), ' in CI = ', num2str(CI_fun(input))])
plotinput2b = input; 
disp(['----------'])      


%Model:
input = nanmean(s0HjSQ,2)*100;
plotinput3a = input; 
disp(['MODEL: Average probability to trade with the home dealer on the platform is (conditional on entry): ', num2str(nanmean(nanmean(s0HjSQ))*100),  ' in CI = ', num2str(CI_fun(input))])

%Data:
data = Y(Y(:,4)==retailer & Y(:,5)==side,:); 
s0HjSQ_data = zeros(T,max(Dmax));
for t=1:T
    for j=1:max(Dmax)          
            N_H = sum(data(data(:,3)==t & data(:,17)==j & data(:,10)==2,26));
            N_T = sum(data(:,3)==t & data(:,17)==j & data(:,10)==2);
            s0HjSQ_data(t,j) = N_H/N_T;
    end 
end  
input = nanmean(s0HjSQ_data,2)*100;
plotinput3b = input; 
disp(['DATA: Average probability to trade with the home dealer on the platform is (conditional on entry): ', num2str(nanmean(input)),  ' in CI = ', num2str(CI_fun(input))])
disp(['----------'])  


%Model
input = nanmean(1-rhojSQ,2)*100;
plotinput4a = input; 
disp(['MODEL: Average probability to enter the platform: ', num2str(nanmean(nanmean(1-rhojSQ))*100), '%,  in CI = ', num2str(CI_fun(input))])

%Data:
data = Y(Y(:,4)==retailer & Y(:,5)==side,:); 
rhojSQ2_data = zeros(T,max(Dmax));
for t=1:T
    for j=1:max(Dmax)          
            N_H = sum(data(:,3)==t & data(:,17)==j & data(:,10)==1);
            N_T = sum(data(:,3)==t & data(:,17)==j);
            rhojSQ2_data(t,j) = N_H/N_T;
    end 
end 
input = nanmean(1-rhojSQ2_data,2)*100;
plotinput4b = input; 
disp(['DATA: Average probability to enter the platform: ', num2str(nanmean(input)), '%,  in CI = ', num2str(CI_fun(input))])
disp(['----------'])  


%Platform cost 
disp(['MODEL: Average platform cost is:', num2str(-1*mean(Y(new_period(:,2),30)))])
disp(['DATA: Average platform cost is 1.1-1.5 bps'])
plotinputcosta = -Y(new_period(:,2),30); 
l1=ceil(size(plotinputcosta)/2);
l2=floor(size(plotinputcosta)/2);
plotinputcostb = [1.5*ones(l1(1),1); 1.1*ones(l2(1),1)];
    
    
%PLOT DISTRIBUTIONS
X=[plotinput1a, plotinput1b, plotinput2a*ones(size(plotinput2b)), plotinput2b, plotinput4a, plotinput4b, plotinput3a, plotinput3b ];
color2 =[ones(size(X,2), 1), ones(size(X,2), 1), ones(size(X,2), 1)];
color1 =[zeros(size(X,2), 1),zeros(size(X,2), 1), zeros(size(X,2), 1)] ; 

figure 
x = 1:size(X,2);
boxplot(rmoutliers(X,'percentiles',[1,99]) , x, 'Symbol','+k', 'Color', 'k', 'orientation', 'vertical', ...
        'Label', {'Pr(H|I)',  '','Pr(H|R)',  '', 'Pr(P|I)','','Pr(H|I,P)' ,''})
ylim([0 110])
h = findobj(gca,'Tag','Box');
for j=1:length(h)
    if rem(j,2)==1
    patch(get(h(j),'XData'),get(h(j),'YData'),color1(j,:),'FaceAlpha',.1);
    else 
    patch(get(h(j),'XData'),get(h(j),'YData'),color2(j,:),'FaceAlpha',.3);
    end
end
ylabel('%')
set(gca,'FontSize',20,'XTickLabelRotation',0)
bp = gca;
set(gcf,'color','w');
temp=['/modelfit0',['_',num2str(side),'_' num2str(retailer),'_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(tradesize), '_', num2str(fej)]]; 
saveas(gcf,fullfile([folderfig,temp,'.png']))


disp(['MARKET CONCENTRATION--------------------------------------------------------------------------'])
disp(['----------'])  

%Model
%platform market share split into home dealer and dealer
for j=1:9
    sNHjSQ_total(:,j) =nansum(sNHjSQ(:,:,j),2);
    s0NHjSQ_total(:,j)=nansum(s0NHjSQ(:,:,j),2);
end 

%HHI on the paltform 
%Model
s0jSQ_total          = s0NHjSQ_total + s0HjSQ;
mass_investors_plat  = nansum(s0jSQ_total,2);
market_share_plat    = (s0jSQ_total)./mass_investors_plat; 
HHI_plat             = HHI_fun(market_share_plat);

input = HHI_plat;
plotinput5a = input; 
disp(['MODEL: Average HHI on the paltform is: ', num2str(nanmean(HHI_plat)), ' in CI = ', num2str(CI_fun(input))])

%Data 
data = Y(Y(:,4)==retailer & Y(:,5)==side,:); 
s0jSQ_total_data = zeros(T,max(Dmax));
for t=1:T
    for j=1:max(Dmax)          
            N_H = sum(data(:,3)==t & data(:,17)==j & data(:,10)==2);
            N_T = sum(data(:,3)==t & data(:,10)==2);
            s0jSQ_total_data(t,j) = N_H/N_T;
    end 
end
mass_investors_plat_data  = nansum(s0jSQ_total_data,2);
market_share_plat_data    = (s0jSQ_total_data)./mass_investors_plat_data; 
HHI_plat_data             = HHI_fun(market_share_plat_data);
input = HHI_plat_data;
plotinput5b = input; 
disp(['DATA: Average HHI on the paltform is: ', num2str(nanmean(input)), ' in CI = ', num2str(CI_fun(input))])
disp(['----------'])  
     
%HHI over-all, including retailers 
%Model
sjSQ_total     = sNHjSQ_total + sHjSQ;
mass_investors = nansum(sjSQ_total + rhojSQ,2);
market_share   = ((1-kappa)*(sjSQ_total + rhojSQ) + kappa)./mass_investors; 
HHI = HHI_fun(market_share);

input = HHI;
plotinput6a = input; 
disp(['MODEL: Average HHI overall is: ', num2str(nanmean(input)), ' in CI = ', num2str(CI_fun(input))])

%Data 
data = Y(Y(:,5)==side,:); 
market_share_data = zeros(T,max(Dmax));
for t=1:T
    for j=1:max(Dmax)          
            N_H = sum(data(:,3)==t & data(:,17)==j);
            N_T = sum(data(:,3)==t);
            market_share_data(t,j) = N_H/N_T;

    end 
end
HHI_data = HHI_fun(market_share_data);

input = HHI_data;
plotinput6b = input; 
disp(['DATA: Average HHI overall is: ', num2str(nanmean(input)), ' in CI = ', num2str(CI_fun(input))])
disp(['----------']) 



%PLOT DISTRIBUTIONS
X = [plotinput6a, plotinput6b, plotinput5a, plotinput5b]; 
color2 =[ones(size(X,2), 1), ones(size(X,2), 1), ones(size(X,2), 1)];
color1 =[zeros(size(X,2), 1),zeros(size(X,2), 1), zeros(size(X,2), 1)] ; 

figure 
x = 1:size(X,2);
boxplot(rmoutliers(X,'percentiles',[1,99]) , x, 'Symbol','+k', 'Color', 'k', 'orientation', 'vertical', ...
        'Label', { 'HHI', '',  'P-HHI' ,''})
ylim([0 1])
h = findobj(gca,'Tag','Box');
for j=1:length(h)
    if rem(j,2)==1
    patch(get(h(j),'XData'),get(h(j),'YData'),color1(j,:),'FaceAlpha',.1);
    else 
    patch(get(h(j),'XData'),get(h(j),'YData'),color2(j,:),'FaceAlpha',.3);
    end
end
set(gca,'FontSize',20,'XTickLabelRotation',0)
bp = gca;
set(gcf,'color','w');
temp=['/modelfitMC',['_',num2str(side),'_' num2str(retailer),'_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(tradesize), '_', num2str(fej)]]; 
saveas(gcf,fullfile([folderfig,temp,'.png']))



disp(['YIELDS-------------------------------------------------------------------------------------'])
disp(['----------'])     
    
%%RETAIL INVESTORS
retailer=1;
for side=sideshow

    %load CF 
    CF=0; 
        %This CF at fixed quotes uses kappa =1 (only retail investors) to
        %compute the expected yield for only retail investors and not a representative investors (this used to be kappa=1, CF=1 before August 2022)
    
    if zerobenefit==0  
        temp=['/counterfac_',[num2str(side), '_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(0),'_' num2str(CF),'_',num2str(flexible), '_', num2str(fej), '_', num2str(loyalty)]];
    else 
       temp=['/counterfac_',[num2str(side), '_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(0),'_' num2str(CF),'_',num2str(flexible), '_', num2str(fej), '_', num2str(loyalty),  '_zerobenefit']];
       homedealer=0;
    end 
    load(fullfile([foldermats,temp,'.mat']))
    Yield_access =Yield21;
 
    %load SQ 
    if zerobenefit==0
        temp=['/counterfac_SQ_',[num2str(side),'_' num2str(retailer), '_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(0), '_', num2str(fej), '_', num2str(loyalty)]]; 
    else
        temp=['/counterfac_SQ_',[num2str(side),'_' num2str(retailer), '_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(0), '_', num2str(fej), '_', num2str(loyalty), '_zerobenefit']]; 
        homedealer=0;
    end
    load(fullfile([foldermats,temp,'.mat']))


    %Reduction in yield when losing platform access
    if side==1
        YieldReduction=-(Yield_access-Yno);
        input=-nanmean(YieldReduction,2);
        disp(['MODEL: Without platform access, retail investors expect a yield that is lower by ', num2str(nanmean(input)), ' bps when buying in CI = ', num2str(CI_fun(input))])
        disp(['DATA: see event study'])
        inputplotYDa=input;
        
        muexp = 1.160;
        stderr = 0.340;      
        inputplotYDb=normrnd(muexp, stderr, [T,1]);
       
    elseif side==-1
        YieldReduction=(Yield_access-Yno);
        input=nanmean(YieldReduction,2);
        disp(['MODEL: Without platform access, retail investors expect a yield that is lower by ', num2str(nanmean(nanmean(YieldReduction))), ' bps when selling in CI = ', num2str(CI_fun(input))])
        disp(['DATA: see event study'])
        
        inputplotYDa=input;
         muexp = 1.160;
         stderr = 0.340;      
         inputplotYDb=normrnd(muexp, stderr, [T,1]);
    end 


end 

%%INSTITUTIONAL INVESTORS
retailer=0;
side=sideshow;

    %load SQ 
    if zerobenefit==0
        temp=['/counterfac_SQ_',[num2str(side),'_' num2str(retailer), '_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(tradesize), '_', num2str(fej), '_', num2str(loyalty)]]; 
    else
        temp=['/counterfac_SQ_',[num2str(side),'_' num2str(retailer), '_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(tradesize), '_', num2str(fej), '_', num2str(loyalty), '_zerobenefit']]; 
        homedealer=0;
    end
    load(fullfile([foldermats,temp,'.mat']))
    Yield_access=YieldSQ;

    %Reduction in yield when losing platform access
    if side==1
        YieldReduction=-(Yield_access-Yno);
        input=nanmean(YieldReduction,2);
        disp(['MODEL: Without platform access, instititional investors expect a yield that is lower by ', ...
        num2str(nanmean(nanmean(YieldReduction))), ' bps when buying  in CI = ', num2str(CI_fun(input))])
        disp(['DATA: see event study'])
    elseif side==-1
        YieldReduction=(Yield_access-Yno);
        input=nanmean(YieldReduction,2);
        disp(['MODEL: Without platform access, instititional investors expect a yield that is lower by ', ...
        num2str(nanmean(nanmean(YieldReduction))), ' bps when selling  in CI = ', num2str(CI_fun(input))])
        disp(['DATA: see event study'])
    end 

    %Platform yield vs. bilateral yield
    %Data:
    data = Y(Y(:,4)==retailer & Y(:,5)==side,:); 
    YB_data = zeros(T,1);
    YE_data = zeros(T,1);
  
    
    for t=1:T
             YB_data(t) = nanmean(data(data(:,3)==t & data(:,10)==1,9));
             YE_data(t) = nanmean(data(data(:,3)==t & data(:,10)==2,9));
    end 

    %Model
    if side==1
        input=nanmean(q_data-YBSQ./rhojSQ, 2);
        disp(['MODEL: The difference between the expected bilateral and platform yield is: ',   num2str(nanmean(input)), ' in CI = ', num2str(CI_fun(input))])
        plotinput7a = input; 
        input=nanmean(YE_data-YB_data,2);
        disp(['DATA: The difference between the expected bilateral and platform yield is: ',   num2str(nanmean(input)), ' in CI = ', num2str(CI_fun(input))])
        plotinput7b = input; 
    else 
         input=nanmean(-q_data+YBSQ./rhojSQ, 2);
         disp(['MODEL: The difference between the expected bilateral and platform yield is: ',   num2str(nanmean(input)), ' in CI = ', num2str(CI_fun(input))])
         plotinput7a = input; 
         input=nanmean(-YE_data+YB_data,2);
         disp(['DATA: The difference between the expected bilateral and platform yield is: ',   num2str(nanmean(input)), ' in CI = ', num2str(CI_fun(input))])
         plotinput7b = input; 

    end


    %Model
    %load SQ for retailers
    if zerobenefit==1
        temp=['/counterfac_SQ_',[num2str(side),'_' num2str(1), '_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(0), '_', num2str(fej), '_', num2str(loyalty), '_zerobenefit']]; 
    else
        temp=['/counterfac_SQ_',[num2str(side),'_' num2str(1), '_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(0), '_', num2str(fej), '_', num2str(loyalty)]]; 
        load(fullfile([foldermats,temp,'.mat']))
    end 
    YieldSQ_R=YieldSQ;
    
    %load SQ for institutionals
    if zerobenefit==1
        temp=['/counterfac_SQ_',[num2str(side),'_' num2str(0), '_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(tradesize), '_', num2str(fej), '_', num2str(loyalty), '_zerobenefit']]; 
    else
        temp=['/counterfac_SQ_',[num2str(side),'_' num2str(0), '_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(tradesize), '_', num2str(fej), '_', num2str(loyalty)]]; 
    end 
        load(fullfile([foldermats,temp,'.mat']))
    
    YieldSQ_I=YieldSQ;
    if side==1
        Yield_gap_SQ    = YieldSQ_I-YieldSQ_R;
    else 
        Yield_gap_SQ    = YieldSQ_R-YieldSQ_I;
    end 
    input=nanmean(Yield_gap_SQ,2);
    plotinput8a= input; 
    disp(['MODEL: The yield gap between R and I in the SQ is: ' , num2str(nanmean(input)), ' bsp in CI = ', num2str(CI_fun(input))])

    %Data:
    data = Y(Y(:,5)==side,:); 
    YR_data = zeros(T,1);
    YI_data = zeros(T,1);
  
    for t=1:T
             YI_data(t) = nanmean(data(data(:,3)==t & data(:,4)==0,9));
             YR_data(t) = nanmean(data(data(:,3)==t & data(:,4)==1,9));
    end 
    
    if side==1
        Yield_gap_data = YI_data - YR_data;
    else 
        Yield_gap_data = YR_data - YI_data;
    end 
     input=nanmean(Yield_gap_data,2);
    plotinput8b= input; 
    disp(['DATA: The yield gap between R and I in the SQ is: ' , num2str(nanmean(input)), ' bsp in CI = ', num2str(CI_fun(input))])

    
%PLOT DISTRIBUTIONS
X = [plotinput8a, plotinput8b,  plotinput7a, plotinput7b, inputplotYDa, inputplotYDb]; %--> vertical
color2 =[ones(size(X,2), 1), ones(size(X,2), 1), ones(size(X,2), 1)];
color1 =[zeros(size(X,2), 1),zeros(size(X,2), 1), zeros(size(X,2), 1)] ; 

figure 
x = 1:size(X,2);
boxplot(rmoutliers(X,'percentiles',[1,99]) , x, 'Symbol','+k', 'Color', 'k', 'orientation', 'vertical', ...
        'Label', {'I-R-Gap', '', 'P-B-Gap' , '', 'Y-Drop', ''})
h = findobj(gca,'Tag','Box');
for j=1:length(h)
    if rem(j,2)==1
    patch(get(h(j),'XData'),get(h(j),'YData'),color1(j,:),'FaceAlpha',.1);
    else 
    patch(get(h(j),'XData'),get(h(j),'YData'),color2(j,:),'FaceAlpha',.3);
    end
end
ylim([-3,10])
ylabel('bps') 
set(gca,'FontSize',20,'XTickLabelRotation',0)
bp = gca;
set(gcf,'color','w');
temp=['/modelfit2',['_',num2str(side),'_' num2str(retailer),'_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(tradesize), '_', num2str(fej)]]; 
saveas(gcf,fullfile([folderfig,temp,'.png']))
    


    
    